clear all
clc

load('Results1')

%% Employment

LLBas                = lLBas;
LLBas(:,2:end,2:end) = repmat(lLBas(:,2:end,1),1,1,T-1).*cumprod(LDBas,3);
LDCfa                = LHCfa .* LDBas;
LLCfa                = lLCfa;
LLCfa(:,2:end,2:end) = repmat(lLCfa(:,2:end,1),1,1,T-1).*cumprod(LDCfa,3);
PopulationB          = squeeze(sum(lLBas(1:M,:,:),2));
TotalEmploymentB     = squeeze(sum(LLBas(1:M,2:end,:),2));
PopulationC          = squeeze(sum(lLCfa(1:M,:,:),2));
TotalEmploymentC     = squeeze(sum(LLCfa(1:M,2:end,:),2));
LoPB                 = TotalEmploymentB./PopulationB;
LoPC                 = TotalEmploymentC./PopulationC;
LoPchan              = (mean(LoPC(:,7:9),2)-mean(LoPB(:,7:9),2))*100;
disp('S.d. of changes in employment-to-pop ratio mentioned in Section 6.1')
disp(compose('%.2f',std(LoPchan)))

%% Income

% Compute things in baseline across three years (06, 07, and 08) to average

Population           = squeeze(sum(lLBas(:,:,:),2));
pop2000              = Population(:,1);
Pagg                 = [ones(I,1),squeeze(prod(PDBas .^ ...
                        repmat(Amat,1,1,T-1),2))];
cumP                 = cumprod(Pagg,2);
labinc2006           = sum(YLBas(:,:,7),2);
price2006            = cumP(:,7);
rlabinc2006          = labinc2006./price2006;
pop2006              = Population(:,7);
incperca2006         = rlabinc2006./pop2006;
incperca2006US       = incperca2006(1:50);
labinc2007           = sum(YLBas(:,:,9),2);
price2007            = cumP(:,8);
rlabinc2007          = labinc2007./price2007;
pop2007              = Population(:,8);
incperca2007         = rlabinc2007./pop2007;
incperca2007US       = incperca2007(1:50);
labinc2008           = sum(YLBas(:,:,9),2);
price2008            = cumP(:,9);
rlabinc2008          = labinc2008./price2008;
pop2008              = Population(:,9);
incperca2008         = rlabinc2008./pop2008;
incperca2008US       = incperca2008(1:50);
incpercaavgUSbas     = (incperca2006US+incperca2007US+incperca2008US)/3;

% Compute things in counterfactual

PDCfa                = PHCfa.*PDBas;
Pagg                 = [ones(I,1),squeeze(prod(PDCfa .^ ...
                        repmat(Amat,1,1,T-1),2))];
cumP                 = cumprod(Pagg,2);
Population           = squeeze(sum(lLCfa(:,:,:),2));
labinc2006           = sum(YLCfa(:,:,7),2);
price2006            = cumP(:,7);
rlabinc2006          = labinc2006./price2006;
pop2006              = Population(:,7);
incperca2006         = rlabinc2006./pop2006;
incperca2006US       = incperca2006(1:50);
labinc2007           = sum(YLCfa(:,:,8),2);
price2007            = cumP(:,8);
rlabinc2007          = labinc2007./price2007;
pop2007              = Population(:,8);
incperca2007         = rlabinc2007./pop2007;
incperca2007US       = incperca2007(1:50);
labinc2008           = sum(YLCfa(:,:,9),2);
price2008            = cumP(:,9);
rlabinc2008          = labinc2008./price2008;
pop2008              = Population(:,9);
incperca2008         = rlabinc2008./pop2008;
incperca2008US       = incperca2008(1:50);
incpercaavgUScfa     = (incperca2006US+incperca2007US+incperca2008US)/3;

% Now compute changes between baseline and counterfactual and display std

logchaincperca       = log(incpercaavgUScfa./incpercaavgUSbas);
sharespopUS          = pop2000(1:50)/sum(pop2000(1:50));
sumchainc            = sharespopUS' * logchaincperca;
difference           = (logchaincperca-sumchainc)*100;
disp('S.d. of changes in income per capita mentioned in Section 6.1')
disp(compose('%.1f',std(difference)))
